/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "uhm/common.hxx"
#include "uhm/const.hxx"

#include "uhm/matrix/uhm/matrix.hxx"
#include "uhm/matrix/uhm/fla.hxx"

namespace uhm {
  void Matrix_FLA_::check_lu_nopiv_1() {

#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix 
    // ----------------------------------------------------------
    if (this->fs) {
      FLA_Copy( ~(this->flat.xt), ~(this->flat.rt) );
      linal::trmm( FLA_LEFT, FLA_UPPER_TRIANGULAR, FLA_NO_TRANSPOSE,
                          FLA_NONUNIT_DIAG, FLA_ONE, 
                          this->hier.ATL, this->hier.rt ); 
    }

    if (this->fs && this->ss) {
      linal::gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE, FLA_ONE,
                          this->hier.ATR, this->hier.xb, 
                          FLA_ONE, this->hier.rt );
    }
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix 
    // ----------------------------------------------------------
    if (this->fs) {
      FLA_Copy( ~(this->flat.xt), ~(this->flat.rt) );
      FLA_Trmm( FLA_LEFT, FLA_UPPER_TRIANGULAR, FLA_NO_TRANSPOSE,
		FLA_NONUNIT_DIAG, FLA_ONE, 
		~(this->flat.ATL), ~(this->flat.rt) ); 
    }

    if (this->fs && this->ss) {
      FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE, FLA_ONE,
		~(this->flat.ATR), ~(this->flat.xb), 
		FLA_ONE, ~(this->flat.rt) );
    }
#endif
  }
   
  // from leaf to root
  void Matrix_FLA_::check_lu_nopiv_2() {
#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Flat-Matrix 
    // ----------------------------------------------------------
    if (this->fs && this->ss) {
      linal::gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE, FLA_ONE,
                          this->hier.ABL, this->hier.rt, 
                          FLA_ONE, this->hier.rb );
    }
    
    if (this->fs) {
      linal::trmm( FLA_LEFT, FLA_LOWER_TRIANGULAR, FLA_NO_TRANSPOSE,
                          FLA_UNIT_DIAG, FLA_ONE, 
                          this->hier.ATL, this->hier.rt ); 
    }
    // rb should be merged for upper hierarchy
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix 
    // ----------------------------------------------------------
    if (this->fs && this->ss) {
      FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE, FLA_ONE,
		~(this->flat.ABL), ~(this->flat.rt), 
		FLA_ONE, ~(this->flat.rb) );
    }
    
    if (this->fs) {
      FLA_Trmm( FLA_LEFT, FLA_LOWER_TRIANGULAR, FLA_NO_TRANSPOSE,
		FLA_UNIT_DIAG, FLA_ONE, 
		~(this->flat.ATL), ~(this->flat.rt) ); 
    }
    // rb should be merged for upper hierarchy
#endif
  }
}
